function [S,detF] = solid_angle(V,F,O,legacy)
  % SOLID_ANGLE Compute the solid angle of a triangle (tetrahedron) described
  % by points V
  %
  % Inputs:
  %  V  #V by dim list of vertex positions
  %  F  #F by #simplex_size list of triangle indices
  %  O  #O by dim list of origin positions {[0,0,0]}
  % Outputs:
  %  S  #O by #F list of solid angles
  %
  warning('Compile gptoolbox/mex/solid_angle.cpp for better performance ...');

  dim = size(V,2);
  assert(size(F,2) == dim);

  if ~exist('legacy','var')
    legacy = false;
  end

  switch dim
  case 3
    if legacy
    %if size(O,1) < 1000
      V = permute(bsxfun(@minus,V,permute(O,[3 2 1])),[1 3 2]);

      % http://en.wikipedia.org/wiki/Solid_angle#Tetrahedron

      % determinant for each triangle with vectors as columns
      detF = ...
        V(F(:,1),:,1).*V(F(:,2),:,2).*V(F(:,3),:,3)+ ...
        V(F(:,2),:,1).*V(F(:,3),:,2).*V(F(:,1),:,3)+ ...
        V(F(:,3),:,1).*V(F(:,1),:,2).*V(F(:,2),:,3)- ...
        V(F(:,3),:,1).*V(F(:,2),:,2).*V(F(:,1),:,3)- ...
        V(F(:,2),:,1).*V(F(:,1),:,2).*V(F(:,3),:,3)- ...
        V(F(:,1),:,1).*V(F(:,3),:,2).*V(F(:,2),:,3);
      % edge lengths numbered same as opposite vertices
      NF = cat(3, ...
        sqrt(sum(V(F(:,1),:,:).^2,3)), ...
        sqrt(sum(V(F(:,2),:,:).^2,3)), ...
        sqrt(sum(V(F(:,3),:,:).^2,3)) ...
        );
      S = 2*atan2( ...
        detF, ...
        prod(NF,3) + ...
        sum(V(F(:,1),:,:).*V(F(:,2),:,:),3).*NF(:,:,3) + ...
        sum(V(F(:,2),:,:).*V(F(:,3),:,:),3).*NF(:,:,1) + ...
        sum(V(F(:,3),:,:).*V(F(:,1),:,:),3).*NF(:,:,2) ...
      )';

    else
      S = zeros(size(O,1),size(F,1));

      %% loop over origins
      %for o = 1:size(O,1)
      %  S(:,o) = solid_angle(V,F,O(o,:),true);
      %end

      V = permute(bsxfun(@minus,V,permute(O,[3 2 1])),[1 3 2]);
      % loop over faces
      for f = 1:size(F,1)
        % determinant for each triangle with vectors as columns
        detF = ...
          V(F(f,1),:,1).*V(F(f,2),:,2).*V(F(f,3),:,3)+ ...
          V(F(f,2),:,1).*V(F(f,3),:,2).*V(F(f,1),:,3)+ ...
          V(F(f,3),:,1).*V(F(f,1),:,2).*V(F(f,2),:,3)- ...
          V(F(f,3),:,1).*V(F(f,2),:,2).*V(F(f,1),:,3)- ...
          V(F(f,2),:,1).*V(F(f,1),:,2).*V(F(f,3),:,3)- ...
          V(F(f,1),:,1).*V(F(f,3),:,2).*V(F(f,2),:,3);
        % edge lengths numbered same as opposite vertices
        NF = cat(3, ...
          sqrt(sum(V(F(f,1),:,:).^2,3)), ...
          sqrt(sum(V(F(f,2),:,:).^2,3)), ...
          sqrt(sum(V(F(f,3),:,:).^2,3)) ...
          );
        % Van Oosterom, A; Strackee, J (1983). "The Solid Angle of a Plane
        % Triangle".
        S(:,f) = 2*atan2( ...
          detF, ...
          prod(NF,3) + ...
          sum(V(F(f,1),:,:).*V(F(f,2),:,:),3).*NF(:,:,3) + ...
          sum(V(F(f,2),:,:).*V(F(f,3),:,:),3).*NF(:,:,1) + ...
          sum(V(F(f,3),:,:).*V(F(f,1),:,:),3).*NF(:,:,2) ...
        );
        %% Eriksson 1990
        %S(:,f) = 2*atan2( ...
        %  detF, ...
        %  1 + ...
        %  sum(V(F(f,1),:,:).*V(F(f,2),:,:),3)+ ...
        %  sum(V(F(f,2),:,:).*V(F(f,3),:,:),3)+ ...
        %  sum(V(F(f,3),:,:).*V(F(f,1),:,:),3)...
        %);
      end
    end
  case 2
    % source and destinations
    VS = V(F(:,1),:);
    VD = V(F(:,2),:);
    % 2d vectors from O to VS and VD
    O2VS = bsxfun(@minus,permute(O(:,1:2),[1 3 2]),permute(VS(:,1:2),[3 1 2]));
    O2VD = bsxfun(@minus,permute(O(:,1:2),[1 3 2]),permute(VD(:,1:2),[3 1 2]));
    % WHY WAS THIS NORMALIZATION HERE?! THIS SHOULD NOT BE NECESSARY
    %%% normalize
    %%O2VS = bsxfun(@rdivide,O2VS,sqrt(sum(O2VS.^2,3)));
    %%O2VD = bsxfun(@rdivide,O2VD,sqrt(sum(O2VD.^2,3)));
    %%% http://stackoverflow.com/questions/2150050/finding-signed-angle-between-vectors
    %%%O2VS
    %%% O2VD
    S = -atan2( ...
      O2VD(:,:,1).*O2VS(:,:,2) - O2VD(:,:,2).*O2VS(:,:,1), ...
      O2VD(:,:,1).*O2VS(:,:,1) + O2VD(:,:,2).*O2VS(:,:,2) );
    %S = atan2( ...
    %  (-1-V(:,1)).*(0-V(:,2)) - (0-V(:,2)).*(1-V(:,1)), ...
    %  (-1-V(:,1)).*(1-V(:,1)) + (0-V(:,2)).*(0-V(:,2)) );
  otherwise
    error(sprintf('Unsupported dimension %d',dim));
  end
end
